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1 Main result 

We consider a discrete-time dynamical system with noisy observations 


x t +i = f(x t ,u t ,9) 

Yt ~p Xt (yt ) 


where x t £ R" denotes the system’s state, u t £ R m is a sequence of inputs to be designed and 0 € R p is a 
vector of unknown parameters that we wish to estimate. Observations are drawn independently from a known 
distribution that is parametrized by the system state x t . We assume that for all x t £ R n the probability 
distribution is absolutely continuous with respect to some measure p and we denote its density with respect 
to fi by Px t (yt)- We further assume that this density is differentiable with respect to the parameter Xt and 
define the Fisher information matrix as 





logp Xt (F t )) (V xt log Ps t (F t )) 



We consider this system over a finite horizon 0 < t < N. Our goal is to design a sequence u that provides 
a maximal amount of information about the unknown parameter vector 6 for 6 in a neighbourhood of some 
nominal value of the parameters 9$. Mathematically, we wish to choose u to maximize a function <P(Xy(9o)) 
of the Fisher information that the joint output Y = (To, • ■ •, Yn) carries about the parameter 9. The function 
(j> is chosen to be a measure of the “largeness” of the positive semidefinite matrix X. Multiple choices for the 
function </> have been proposed [T|; here we use 4>{X) = tr(X) (often called “T-optimal design”). In general 
this problem is nonconvex as a function of it, but due to the fact that the trace is linear our objective function 
is additive so a global solution can be found using dynamic programming. 

The following result allows us to compute the information contained in the observed data. A proof of 
this proposition is given in Section [3] 


Proposition 1. Suppose that for all x t £ 
there exists a /. i-integrable function q with 
Fisher information with respect to the parameter 9 can he computed as 


the density p Xt (y t ) is differentiable with respect to x t 
< q(yt ) for all yt £ R. If f is C 1 in Xt and 9, 


dpx t (yt) 


dx t 


and that 
then the 


N 

T y {9) = J2^ sX t) T X Y fftt)(V e x t ) (2) 

t =0 

where VeXt denotes the Jacobian of Xt with respect to 6. 
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Thus, the T-optimal design criterion, the objective function is given by 


N 

tr(2V(0)) = ^tr ((V e x t ) T l Yt (x t )(V 9 x t )y (3) 

t=0 

Applying the chain rule to (JTJ) , we get a dynamical system 

Ygx t+1 = Vgf(x t ,u t ,e) + \7 x f(x t ,u t ,9) \7 e x t (4) 

describing the time evolution of the sensitivities Y gx t - The dynamics ([!]) and Q together with the cost 
function ©> define a discrete-time finite-horizon optimal control problem that can be solved via dynamic 
programming [2]. In particular, we define the sequence of value functions J k via 

■Xv(:EiV, VflTlv) = tr (^CVgXN) T Iy N {XN){'VeXN) S j 

Jk(x k ,Y g x k ) = maxjtr (( VgX k ) T l Yk {x k )(Vgx k )) + J k+1 ^f(x k ,u k ,9),'Vgf(x k ,u k ,6) + V x f(x k ,u k ,9) \7 e x k 

If the control policy u* k = fJ-%{x k , Vgx k ) maximizes the right hand side of (JT|) then u* is globally optimal. 

Related approaches to the optimal experiment design problem appear in d] and [4] for continuous¬ 
time dynamical systems with Gaussian noise. However, a different objective function <f> is used and these 
approaches requires appending a nonlinear matrix differential equation for the dispersion (the inverse of 
the Fisher information) to the system state in addition to equation Q. By choosing the T-optimal design 
criterion 4>{I) = tr(Z), we are able to avoid adding an equation for the dispersion to the system state, 
allowing us to efficiently solve problems of larger dimension. 


2 Example problem 

We consider a population of fruit flies, whose dynamics are modelled using the discrete logistic equation 

x t +i = x t + rx t (K - x t ). 

We want to estimate the reproduction rate r along with the carrying capacity K. To generate data from 
which to estimate the model parameters, we place a sequence of traps into the fly cage, each capturing a 
fraction ut of the current fly population. By measuring the number of fruit flies caught in the trap, we wish 
to infer the model parameters. The optimization problem thus consists of choosing the size of the traps (and 
hence the proportion of flies trapped) at each sampling interval. 

This leads to a model for the population dynamics together with the number of fruit flies trapped Y t 

x 0 = K 

xt+i = x t (l - ut) + rx t ( 1 - u t )(K - Xt (1 - u t )) (5) 

Y t ~ Poisson(xtUt). 

For this problem, we approximate the functions J k by evaluation on a grid of size 100 x 100 x 100. 
We optimize about the nominal parameter values ro = 5 x 10~ 4 and Kq = 1000. This is implemented in 
MATLAB using the dynamic programming routine introduced in (5). The optimal inputs are computed in 
33.01 seconds and are shown in Figure [Ta| The corresponding state trajectory is shown in Figure [Tb] 

We see that the optimal observation scheme is to first capture a large fraction u\ = 0.9795 of the flies 
allowing us to get a reliable estimate for the carrying capacity K. After this, we capture a fraction ut ~ 0.32 
of the flies, just enough to keep the population constant. This provides maximal sensitivity to the growth 
rate r in a neighbourhood of ro- Indeed, if r > ro we will see the the population of flies grow over time, 
whereas if r < ro the population will shrink toward zero. 


2 



(a) Optimal input trajectory for |5| computed 
using dynamic programming 



(b) State trajectory corresponding to the input 
given in Figure 


la 


Figure 1: Numerically-computed solution to the optimal experiment design problem 


3 Proof of Proposition [T] 

Proof. First, note that the hypotheses of this proposition provide sufficient regularity to exchange the order 
of differentiation with respect to 9 j and integration with respect to y t . Therefore for all i = 1,... ,p and 
t = 0,...,N 


Eg 


S log p e (Y t )' 


09 t 


Slog Pe(yt) 

09 t 


Pe(yt)dp{y t ) = 


dpejyt) 

dOi 


dp(y t ) 


d f d 

= wJ M9,)My,) = W, l = 0 


Now for all i,j = 1,... ,p we can compute the th entry of Ig(Y) as 
~0 log p e {Y) Slog pg(Y)~ 


1-y(0)i,j = E 0 

= Eg 


00 j 


09 , 


0logpg(Y t ) Wy, 0logpe{Y s ) \ 

^ w J 


K t =0 


dO 


N 


= E E * 

t =o 
N 

= E E » 

t =0 
N 

= E E * 

t=o 
N 

= E 

t= o 
N 

= E 

t=o 


Slog pe(Y t ) Slog pe(Y t ) 


00i 


09j 
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N N 

E E Ee 

£=0 s=0,s^£ 


Slogpe(Yt) 

09i 


Eg 


Slog pe(Y s ) 

00i 


<9 log p Xt {e)(Yt) 0logp Xtie) (Y t ) 


06i 


00 3 


, Vx t log ^> v Xt log (y t )) 


Sxt' 


'09. 


Eg 


S0, 


(v* t logp Xt (y t )) (v xt log p Xt (Y t )^j 


0x t 

OK, 


0x± . . Sx^ 
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So 

N 

Zy(0) = ^(Vfl* t ) T Jy t (* t )(V 9 * t ). 

t =0 

□ 
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